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Abstract 


A  theoretical  model  that  addresses  thermochemical  erosion  in  gun  tubes  is  presented.  The 
model  incorporates  two  interior  ballistics  codes — XKTC  and  IBHVG2 — and  the  thermochemical 
code  BLAKE  to  provide  the  necessary  state  variables  and  bulk  species  concentrations  of  the  core 
flow  as  input.  This  erosion  model  utilizes  a  Crank-Nicolson  integration  scheme,  with  dynamic 
gridding  capability  to  account  for  material  ablation,  as  well  as  die  addition  of  energy  sources  and 
heat  transfer  augmentation  due  to  surface  deviations.  A  mass  transport  scheme,  utilizing  the 
Lennaid  Jones  6-12  diffusion  model,  enables  individual  species  to  be  transported  to  the  surface 
from  the  core  flow.  Also  fully  coupled  is  a  separate  thermochemical  routine  which  incorporates 
the  NASA  Lewis  database.  The  code  is  written  modularly,  enabling  the  inclusion  and 
modification  of  existing  submodules.  &osion  results  comparing  a  fielded  kinetic  energy  tank 
round  and  a  candidate  next-generation  tank  round  are  presented.  The  thomochemical  effects  at 
the  surface  are  also  shown  and  discussed. 
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1.  Introduction 


Gun  tubes  are  typically  either  erosion  or  fatigue  limited  in  the  number  of  rounds  that  can  be 
effectively  and  safely  fired  from  them.  Direct-fire  weapons  are  typically  erosion-life  limited  due  to 
the  high-performance  requirements  placed  upon  them.  Ammunition  for  the  M256  tank  gun  is 
designed  to  maximize  the  kinetic  energy  of  the  projectile  at  launch.  New  performance  requirements 
have  reintroduced  gun  bore  erosion  as  a  significant  design  issue  because  erosion  will  limit  the 
number  of  rounds  that  can  be  effectively  and  safely  fired  over  the  life  of  the  gun  tube.  Previous 
“solutions”  had  not  identified  the  fundamental  cause  of  the  erosion,  and  some  discrq)ancies  between 
the  flame  temperature  correlations  [1, 2, 3]  and  the  erosivity  were  never  resolved.  Attempts  to  model 
erosion  using  first  principles  have  been  and  are  currently  being  made  [4, 5, 6],  although  it  is  believed 
that  significant  additional  work  is  still  required. 

A  complete  description  of  the  erosion  process  would  include  a  variety  of  mechanisms  produced 
by  thermal,  chemical,  and  mechanical  effects.  Generally,  thermal  erosion  is  driven  by  high 
convective  heating  produced  by  the  propellant  gases  that  heat  the  in-bore  surface  of  the  gun  tube  and 
cause  it  to  melt.  Chemical  erosion  is  caused  by  surface  chemical  reactions  produced  by  the 
interaction  of  the  gun  tube  gases  with  the  in-bore  surface  of  the  tube.  These  reactions  may  arode  the 
in-bore  surface  directly  through  pyrolysis  or  produce  additional  heat  to  augment  the  convective 
heating,  thus  inducing  melting.  The  thermal  and  chemical  erosion  mechanisms  are  common  and 
generally  interrelated,  while  erosion  due  to  mechanical  effects  (projectile  passage)  is  more  isolated 
and  somewhat  independent,  bi  this  current  work,  we  have  focused  on  the  coupled  thermal  and 
chemical  erosion  mechanisms,  while  not  excluding  the  possible  addition  of  the  mechanical 
contribution  in  the  future. 

2.  Model  Description 


The  model  description  is  composed  of  three  fully  coupled  portions  consisting  of  thermal  ablation 
and  heat  transfer/conduction,  mass  transport,  and  the  imbedded  thermochemistry,  as  well  as  the 
intercoupling  of  fliese  portions. 
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Following  a  modular  engineering  approach  resulted  in  the  integration  of  some  of  the  U.S.  Anny 
Research  Laboratory’s  (ARL)  well-developed  tools  as  well  as  the  development  of  new  modules  to 
suppfy  additional  physics — namely,  an  interior  ballistics  code  to  provide  the  necessary  state  variables 
and  gas  velocities  XKTC  [7],  and  a  chemical  equilibrium  code  BLAKE  [8]  coupled  to  an  interior 
ballistic  code,  IBHVG2  [9]  (thus  creating  the  IBBLAKE  [10]  code,  which  provides  ten^oral  species 
concentrations),  as  well  as  the  gun  tube  heat  transfer-conduction  code  XBR-2D  [11, 12, 13]. 

The  current  modeling  effort  considers  a  melt-w^  model  similar  to  that  proposed  by  Caveny  [14], 
with  the  addition  of  surface  chemistry.  Conceptually,  as  shown  in  Figure  1,  the  surface  heats  from 
convection  until  the  chemical  activation  energy  is  overcome;  at  this  point  a  surface  reaction  occurs, 
releasing  additional  energy  into  the  system  and  producing  various  products.  The  reaction  products 
can  either  remain,  as  some  solids,  or  be  removed  from  the  area  as  liquids  or  gases,  thus  resulting  in 
ablation  or  pyrolysis,  respectively.  As  the  surface  regresses,  the  solids  are  refreshed  when  the  virgin 
substrate  material  is  exposed.  The  following  assunptions  have  been  made  in  the  model: 

(1)  One-dimensional  (1-D)  heat  conduction  is  considered. 

(2)  No  subsurface  diffusion  of  species  or  reactions  is  considered. 

(3)  All  surface  liquids  are  removed. 

(4)  All  surface  gas  products  are  removed. 

(5)  There  is  no  feedback  to  the  interior  ballistics. 

(6)  Chemical  energy  released  is  treated  as  a  source  term. 

(7)  Species  are  chemically  frozen  from  core  flow  to  wall. 
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Gas  Phase  Solid  Phase 


Figure  1.  Conceptual  Erosion  Model. 


3.  Ablation-Conduction  Model  and  Computational 
Approach 

3.1  Analytical  Description  of  the  Ablation-Conduction  Model  [15].  The  in-depth 
temperature  response  of  the  unablated  (solid)  material  is  modeled  using  1-D  heat  conduction 
equation  shown  as  follows: 


,  ar  ^  1  8 

P  dt  r^dr 


(1) 


By  setting  P  =  0  or  P  =  1,  the  planar  or  axisymmetric  form  of  the  governing  equation  can  be 
obtained.  In  this  form  of  die  equation,  the  relevant  material  properties  (density  [p],  specific  heat 
[Cp],  and  conductivity  [A:]),  die  conductivity  may  vary,  with  all  the  others  remaining  continuous.  The 
feature  of  variable  conductivity  was  not  utilized  in  the  results  presented  here. 
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At  the  internal  surface  of  the  gun  tube,  two  separate  boundary  conditions  were  applied  depending 
on  whether  melting  of  the  surface  material  was  occurring.  When  the  surface  material  was  below  flie 
melt  temperature,  a  convective  heat  transfer  boundary  condition  was  applied, 

^  ~  Source,  (2) 

with  a  source  term  due  to  the  additional  energy  released  from  surface  reactions.  Using  the  state 
variables  and  gas  velocity  from  the  XKTC  [7]  code,  the  convective  heat  transfer  coefficient,  h,  is 
determined  using  the  correlation  of  Stratford  and  Beavers  [16].  Once  the  surface  of  the  gun  tube 
reaches  the  melt  temperature,  melting  of  the  surface  material  is  assumed  to  occur.  When  fuUy 
molten  (complete  change  of  phase  from  solid  to  liquid),  the  liquid  metal  is  immediately  removed  or 
ablated  by  the  shearing  action  of  the  gun  gases.  Because  surface  material  is  being  removed  during 
the  ablation  process,  the  surface  location  (or  the  regression  rate)  becomes  an  additional  variable. 
During  the  ablation  process,  two  boundary  conditions  are  applied  as  follows: 

T  =  T 

wall  ^  melt 

The  first  equality  of  the  melting  boundary  condition  simply  states  that  the  temperature  at  the 
interface,  while  the  phase  change  is  occurring,  is  equal  to  the  melt  temperature.  The  melt 
temperature  is  assumed  to  be  a  known  material  property.  The  second  equality  of  the  boundary 
condition  is  obtained  from  an  energy  balance  at  the  melt  surface.  In  addition  to  the  heat  transfer  due 
to  conduction  and  convection,  the  additional  energy  required  due  to  the  change  of  phase  from  solid 
to  liquid  Oatent  heat  of  melting,  L)  is  also  included.  This  boundary  condition  allows  the  regression 
rate  of  the  solid  material  to  be  computed,  given  that  the  latent  heat  of  melting,  L,  is  a  known  material 
property. 


To  provide  closure  for  the  in-deplh  temperature  response  of  the  gun  tube,  a  convective  boundary 
condition  is  applied  to  the  outer  surface  of  the  gun  tube. 
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T.)  -  -tf . 


(3) 


However,  because  the  erosion  process  occurs  during  the  first  several  milliseconds  of  the  firing 
process,  the  heat  can  only  penetrate  a  fiuction  of  the  distance  from  the  inner  surface  of  the  gun  tube 
to  the  outer  surface  of  the  gun  tube.  In  tiiis  case,  it  may  not  be  necessary  to  model  the  radial 
temperature  response  of  the  entire  gun  tube.  By  locating  the  outer  edge  of  the  computational  domain 
so  that  the  temperature  response  is  unchanged  during  the  heating  process,  a  smaller  computational 
domain  may  be  analyzed.  In  this  case,  one  of  two  boundary  conditions  may  be  applied — a  constant 
temperature  boundary  condition  or  an  adiabatic  wall  boundary  condition. 

The  appropriate  depth  of  the  computational  domain  can  be  estimated  using  a  depth  of  penetration 
analysis  [17, 18].  The  depth  of  poietration,  5,  gives  the  approximate  distance  that  the  heat  would 
pmetrate  in  a  given  time  t  and  is  a  function  of  the  thermal  diffusivity  of  the  material,  a. 

6  =  y/l2  a  t.  (4) 

Using  a  factor  of  safety  of  two  or  diree  5  will  place  die  outer  boundary  of  the  computational  domain 
far  enough  from  the  bore  surface  of  the  gun  tube  so  that  the  temperature  at  the  outer  boundary 
remains  unaffected  by  the  heating  of  the  gun  tube  during  firing.  This  analysis  enables  either  single¬ 
firing  or  rapid,  continuous-fire  erosion  prediction. 

The  governing  equations  and  boundary  conditions  are  solved  using  a  Crank-Nicolson  finite- 
difference  technique.  Prior  to  the  onset  of  melting,  the  governing  equations  and  boundary  ccmditions 
are  linear,  and  solutions  are  obtained  in  a  direct  (noniterative)  fashion.  During  the  melting  process, 
the  equations  become  nonlinear,  since  the  dimensions  of  die  computational  domain  are  coupled  with 
the  regression  rate.  An  iterative  approach  is  utilized  during  melting  to  appropriately  address  the 
nonlinearity. 
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Because  the  boundary  of  the  computational  domain  moves  during  the  erosion  event,  a 
transformed  version  of  the  governing  equation  is  employed.  This  allows  the  equations  to  be  solved 
in  a  fixed  computational  space  even  though  the  physical  boundary  is  moving.  A  gPinpraliyftH 
transformation  between  the  computational  coordinate,  5,  and  the  physical  coordinate,  r,  is  utilized. 
The  transformed  equations  are  shown  as  follows: 


as 


In  this  form,  the  nonlinear  nature  of  the  governing  equation  produced  by  the  moving  boundary  is 
evident  because  the  metric  terms,  and  are  not  constant  and  are  dependent  on  die  erosion  rate 
when  the  grid  is  moving. 

3J2  Comparison  of  Ablation-Conduction  Model  With  Semianalytical  Approach  of  Landau 
(Without  Source  Term).  Over  the  past  several  decades,  numerous  studies  examining  the  phase 
change  process  have  been  made.  Results  fi-om  two  of  these  studies  have  been  compared  with  the 
currrat  numerical  method  without  the  inclusion  of  surface  chemistry  for  simplicity.  T  [19]  hag 
made  1-D  time-dependent  numerical  predictions  of  a  melting  solid  using  a  technique  similar  to  that 
utilized  here.  Landau  considered  the  case  of  a  semi-infinite  solid  subject  to  a  constant  heat  flux. 
Goodman  [20]  considered  sumlar  problems  using  a  heat-balance  integral  approach  that  utilizes  an 
assumed  form  of  the  temperature  profile  to  analytically  determine  the  heat  conduction  and  ablation 
process.  While  exact  solutions  are  not  typicaUy  obtained,  the  results  are  reasonably  accurate  and  are 
of  a  simple  form.  One  of  the  cases  addressed  by  Goodman  was  the  melting  of  a  semi-infinite  solid 
under  constant  heat  flux.  In  both  cases,  an  ablative  boundary  condition  was  utilized;  th^t  is,  the 
liquid  phase  was  immediately  removed  following  melting. 
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The  problem  of  the  heating  on  a  semi-iniinite  solid  subject  to  constant  heat  flux  was  addressed 
using  the  current  numerical  method.  The  solid  initially  had  a  uniform  temperature  distribution  of 
Tg  and  a  specified  melt  temperature,  The  temperature  change  is  normalized  by  -  T^,  while 
the  in-wall  depth  is  normalized  by  Figure  2  shows  the  in-depth  temperature  profile  at 

several  increments  in  time  during  the  melting  process.  The  results  were  obtained  using  material 
properties  corresponding  to  gun  steel.  Li  addition  to  the  numerical  results,  the  analytical  results 
obtained  using  the  heat-balance  integral  approach  of  Ckxxhnan  are  also  shown  at  the  onset  of  melting 
and  dining  steady-state  melting.  The  numerical  solution  and  the  analytical  results  are  in  excellent 
agreemrat  at  the  onset  of  melting.  (Though  not  shown,  both  results  are  also  in  excellent  agreement 
with  the  exact  analytical  solution.)  As  the  melting  progresses,  the  temperature  gradient  at  the  heated 
interface  is  reduced  due  to  the  additional  energy  required  to  melt  the  solid  (latent  heat  of  melting) 
as  shown  in  equation  4  repeated  here  for  clarity  without  the  source  term 


pL|  =  *(r, 


gas 


After  approximately  eight  times  the  time  required  to  reach  the  melt  temperature,  the  numerical 
results  show  that  temperature  profile  has  nearly  reached  the  steady-state  temperature  profile.  The 
differences  between  the  analytical  and  numerical  steady-state  temperature  profile  are  less  than  2% 
of  the  difference  between  the  melt  temperature  and  the  initial  temperature.  Exact  agreement  between 
the  two  results  is  not  expected  because  the  assumed  form  of  the  temperature  profile  in  the 
heat-balance  integral  approach  is  not  an  exact  solution. 

A  comparison  of  the  predicted  ablation  rate  obtained  using  the  current  technique  with  the  results 
of  Landau  is  shown  in  Figure  3.  The  current  results  are  in  excellent  agreement  with  the  results  of 
Landau  for  three  different  sets  of  material  properties  including  properties  close  to  gun  steel,  m  =  2.26 
in  Figure  3.  (Differences  between  the  results  of  Landau  and  the  current  result  may  be  more  related 
to  problems  in  extracting  the  data  from  the  original  published  graphs  of  Landau  than  to  numerical 
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Figure  2.  In-Depth  Temperature  Profile  During  Ablation  Process  for  a  Semi-Infinite  Slab 
Subject  to  Constant  Heat  Hux. 


Figure  3.  Fractional  Ablation  Rate  vs.  Time  for  a  Semi-Infinite  Slab  Subject  to  Constant 
Heat  Flux  (m  «  2.26  for  Gun  Steel). 
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accuracy.)  The  steady-state  ablation  rate  for  the  current  result  is  also  in  good  agreement  with  the 
exact  analytical  results  presented  by  Landau. 

Up  to  this  point  in  the  report,  the  convective  energy  transport,  heat  conduction,  and  surface  melt 
wipe  model  have  been  developed.  The  following  section  4  deals  with  die  transport  of  the  mass  from 
the  core  flow  (center  of  the  bore/chamber)  to  the  wall  through  the  turbulent  boundary  layer.  The 
species  tran^orted  to  the  surface  will  then  be  reacted  with  the  surface  through  a  procedure  outlined 
in  section  5. 

4.  Mass  Transport  and  Multicomponent  Diffusion 

The  mass  transport  is  modeled  using  mass  transport  coefficients  not  unlike  the  heat  transport 
coefficient.  However,  because  each  species  may  diffuse  differently  due  to  its  molecular  weight,  the 
proportions  at  the  wall  may  be  different  than  in  the  core  flow.  This  requires  each  species  to  be 
transported  individually  to  the  wall  through  the  mixture.  Mass  transport  coefficients  are  represented 
in  Sherwood  number  correlations  (Sh  =  hJUDj^^,  as  h^,  where  L  is  a  length  parameter,  and  the 
diffusion  coefficient,  Aib,  must  also  be  determined. 

A  general  outline  of  the  mass  transport  procedure  is  as  follows: 

(1)  Determine  the  binaiy  diffusion  coefficient  for  each  q)ecies  using  the  Lennard  Jones  potential 
parameters. 

(2)  Determine  the  mixture  diffusion  coefficient  and  viscosity  for  each  specie  through  Wilke’s 
mixing  rule. 

(3)  Determine  the  mass  transport  coefficient  from  the  Sherwood  correlations. 

(4)  Solve  the  boundary  condition  for  each  species,  similarly  to  the  energy  transport. 
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4.1  Section  1.  A  reasonable  representation  of  binaiy  diffusion  of  species  “1”  into  species  “2” 
is  provided  by  the  following  relationship  using  the  Lennard  Jones  6-12  potential  function 
paramet^  [21]: 


D 


12 


^  0.0026280  (Mj  + 


Pa 


2  na.!)*/-,*) 
12  “i2  ynf 


(5) 


where  Mj,  are  the  molecular  weights  of  the  species,  Tis  temperature,  P  is  pressure,  and,  0,2 ’s  are 
the  collision  diameters  for  low-energy  collisions.  The  collisional  cross-section  integral  quantities, 
Q*i2’s,  are  obtained  through  table  interpolation.  The  potential  function  for  a  particular  species ‘T’ 
governing  the  Leimard  Jones  model  [21]  is  presented  as 


Y.(r)  =  4e. 


12 

6 

1  ’’ij 

(6) 


where  e,  is  the  depth  of  the  potential  weU.  The  potential  weU  depths,  e/s,  and  the  collisional 
diameters,  o/s,  are  provided  in  tabulated  form  for  various  species. 


The  following  nondimensional  temperature  Tij*  is  used  as  a  tabulation  parameter  for  the 
Lennard-Jones  collisional  cross  section,  Q*i2»  integrals,  where  k  is  the  Boltzmans  constant. 


Combinatory  analysis  [21]  is  needed  as  shown  to  correctly  apply  the  tabulated  single  species  values 
in  equation  6  and  nondimensionalize  the  temperature  in  equation  7. 


e 


12 


'12 


(Oi  -*•  O2) 
2 
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4.2  Section  2.  The  computed  binaiy  diffusion  provides  the  basis  for  the  multicomponent 
diffusion  coefficient  Each  binary  diffusion  possibility  is  used  and  weighted  vs.  all  other  possibilities 
in  the  following  mixture  coefficient  combinatory  function  [22]. 


1-X, 


where  X;  is  the  mass  fraction  of  series.  This  function  enables  the  calculation  of  the  diffusion 
coefficient  for  a  particular  specie ‘T’  into  a  mixture  “m”  of  many  species.  Equilibrium  calculations 
for  all  service  propellants  were  performed,  and  all  principle  products  species  o’s  and  k/€'s  along  with 
the  collisional  cross  sections,  QfT*),  were  placed  into  data  statements  in  the  diffusion  routine.  Thus, 
tracking  of  an  individual  species  contribution  to  the  erosion  process  is  possible.  The  aforementioned 
methodology  was  independently  checked  against  known  diffusivities  with  good  agreement. 

Utilizing  the  collisional  cross  sections,  Q^’s,  for  viscocity  and  collisional  diameters,  a’s,  and  the 
molecular  weight,  the  following  relationship  derived  from  kinetic  theory  [21]  is  utilized  to  determine 
the  mixture  viscocity  using  Wilke’s  rule.  The  mixture  viscosity  is  used  in  developing  the  Reynolds 
number  and  Schmidt  number.  Sc  -  p/pl>^,  where  p  is  the  density. 


P 


=  2.6693  X  10-5 


43  Section  3.  The  development  of  the  mass  transport  coefficients  may  now  proceed  with  the 
previously  developed  computational  methodology  for  both  the  diffusion  coefficients  and  viscosities 
of  the  species  in  the  mixture. 
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Mass  transport  to  the  surface  from  the  core  flow  is  provided  through  a  species  concentration 
potential  (i>jcorefltw  ~  ^iwau)  aiid  a  mass  transport  coefficient,  for  each  species.  The  species 
concentration  of  die  core  flow  is  provided  by  IBBLAKE  [10]  while  the  species  concentration  at  the 
surface  “wall”  is  what  is  being  sought. 

Development  of  appropriate  Sherwood  number  correlations  requires  boundary  layer  analysis  and 
experimental  data.  Mass  transport  for  a  transient  boundary  layer  development  as  shown  in  Figure  4 
has  been  derived  [23]  and  will  be  presented  here  for  a  periodic  transitional  sublayer  boundary  layer 
model  that  assumes  infrequent  penetration  to  the  surface  (thin  wall  layer)  by  eddies  entering  the 
boundary  layer.  This  approach  enables  predictions  of  transfer  rates  matching  expoimental  data  over 
a  wide  range  of  Schmidt  numbers.  The  thin  laminar  wall  layer  beneath  the  boundary  layer  is 
poiodically  replenished  widi  species  from  the  sublayer,  however  at  a  much  slowra*  rate  than  the  bulk 
wall  region.  The  importance  of  the  thin  wall  layer  determines  the  treatment  of  the  model. 


Sublayer  Uavalopmant  phase 


Figure  4.  Idealized  Boundary  Layer  Development  [23]. 

At  moderate  Schmidt  numbers  (10  <  Sc  <  1,000),  the  thickness  of  the  boundary  layer  is  much 
greater  than  the  thickness  of  the  viscous  sublayer.  This  enables  the  accumulation  in  the  wall  layer 
to  be  neglected.  Ruckenstein  [24]  assumed  the  following  linear  pressure  gradient  model  to  govern 
the  mass  transport  to  the  surface  through  the  thin  wall  layer 
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where  k  is  the  conductivity  and  dj  is  the  diickness  of  the  thin  wall  layer.  Utilizing  the  momentum 
integral  method  [23],  the  following  Sherwood  number  has  been  derived: 


Sh 


9.  i  (  .J. 

0.0097jge^°5c^  j  1.10  +  0.445c  ^ 

1  {  -1 
1  +  0.064SC  2  ^  1.10  +  0.44SC  ^ 


-  0.70SC  ^ 


-  0.70SC  « 


For  higher  Schmidt  numbers  (5c  >  1,000)  the  thickness  of  the  boundary  layer  becomes  of  the 
order  of  magnitude  of  the  thin  wall  layer.  The  mass  transport  for  this  range  is  governed  by  the 
penetration  of  eddies  into  the  boundary  layer,  in  this  case  the  wall  layer  [23].  This  description, 
assuming  characteristic  distances  between  eddy  dismptions  is  larger  in  the  wall  layer  than  it  is  in  the 
boundary  layer  and  Blasius  expression  for  the  Motion  factor,  leads  to  the  following  Sherwood 
number 

_9_ 

Sh  =  0.0102 /?c  Sc 

These  expressions  for  die  Showood  number  have  bera  compared  [23]  with  much  experimental 
data  and  agree  very  well  within  the  Schmidt  number  regions  specified. 

AA  Section  4.  The  Sherwood  number  values  and  appropriate  diffusivities  enable  the  calculation 
of  the  mass  transport  coefficient  that  closes  the  boundary  condition.  The  mass  transport  coefficirats, 
A„’s,  are  daived  from  Sherwood  number  correlations,  Sh  =  hJJD^,  and  integrated  over  space  and 
time,  as  shown,  to  provide  the  temporally  varying  ^ecies  at  the  surface. 

Mass,  =  I dAdt. 
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Although  the  gaseous  species  are  assumed  to  not  penetrate  the  surface  at  this  time,  the  diffusion 
module  is  general  enough  to  readily  incorporate  this  possibility  in  the  future. 

At  this  point,  both  the  convective  energy  transport  and  mass  transport  have  been  presented. 
Thus,  the  state  propaties  as  well  as  surface  concentrations  of  specific  species  have  been  described. 
At  this  point,  some  treatment  of  the  solid-gaseous  interaction  can  be  made.  Due  to  the  lack  of 
understanding  of  bofli  the  specific  surface  reactions  and  rates  associated  with  these  reactions  under 
balhstic  conditions,  infimte  rate  kinetics  have  been  applied  in  this  model  as  a  first  approximation  to 
the  surface  chemistry.  A  brief  description  of  the  energy  minimization  methodology  for  eijuilibrium 
reactions  follows. 


5.  Equilibrium  Kinetics 


Equilibrium  chemical  processes  are  considered  to  dominate  whenever  the  characteristic  time  for 
a  “fluid  element”  (discritized  element  in  finite  difference  model  of  a  continuum  region  of  a  fluid 
dynaimc  event)  to  traverse  the  “flow  field  of  interest”  (a  fluid  dynamic  event)  is  much  longer  than 
the  characteristic  time  for  chemical  reactions  to  approach  equilibrium.  As  the  pressure  and 
temperature  increase,  the  molecular  collision  frequency  and  energy  per  collision  increase,  which 
leads  to  smaller  characteristic  chemical  times.  Thus,  chemical  processes  approach  equilibrium 
(i.e.,  occur  at  almost  an  infinite  rate). 


Chemical  equilibrium  is  usually  described  by  either  of  two  equivalent  formulations:  equilibrium 
constants  or  nunimization  of  fiiee  aiergy.  Several  disadvantages  of  the  equilibrium  constant  method 
have  been  noted  [25],  and  most  researchers  prefer  the  fi:ee-energy  minimization  formulation.  The 
condition  for  equilibrium  may  be  stated  in  terms  of  any  of  several  thermodynamic  functions  such  as 
the  minimization  of  the  Gibbs  fi«e  energy  or  Helmholtz  free  «iergy  or  the  maximization  of  entropy. 
For  a  mixture  of  N  species  (e.g.,  atoms  or  molecules)  with  the  number  of  moles  of  species  denoted 
Hf,  the  Gibbs  energy  per  mole  of  mixture  is  given  in  terms  of  the  Gibbs  fi^e  energy  of  the  individual 
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species,  the  internal  energy,  e,  the  temperature,  T,  the  entropy,  s,  the  pressure,  p,  and  the  specific 
volume,  V, 

G  =  n.g.  =  e  -  Ts  +  pV. 

The  equilibrium  method  employed  in  the  present  study  is  based  on  the  fact  that  at  equilibrium 
the  total  Gibbs  energy  of  die  system  attains  a  minimum  value.  The  problem  is  to  find  the  set  of  n/s 
that  minimizes  the  Gibbs  free  energy  for  a  specified  energy  and  volume  (e,V),  subject  to  the 
constraints  of  material  balances.  The  standard  solution  to  this  type  of  problem  is  based  on  the 
metiiod  of  Lagrange’s  undetermined  multipliers.  First,  we  must  recognize  that  the  total  number  of 
atoms  of  each  element  in  the  syston  is  constant.  Denoting  a  particular  atomic  species  with  subscript 
k,  then  A*  is  the  total  number  of  atomic  masses  of  the  Ath  element  in  the  system,  as  determined  by 
the  initial  constitution  of  the  system.  Denoting  the  number  of  atoms  of  the  ifcth  element  present  in 
each  molecule  of  chemical  species  i  by  % ,  then  the  material  balance  on  each  element  k  may  be 
written  (Af  used  here  is  the  number  of  elements), 

h  (E?1i  =  0  (^  =  1,2,...,M)  (8) 

introducing  Lagrange  multipliers,  X* ,  for  each  element  Then  a  new  function,  F,  is  formed  by 
addition  of  the  last  equation  to  G^.  The  frmction  F  is  idmtical  to  since  the  summation  term 
is  zero.  However,  dFldrij  and  dG^drii  are  different  since  F  incorporates  the  constraints  of  the 
matoial  balances.  The  minimum  of  both  F  and  G^  occurs  when  these  partial  derivatives  are  zero. 


F 


e.V,  n, 


dG. 


total 


8n, 


-Er.. 


=  0 


e.V.n, 


UorF^). 
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This  equation  can  be  rewritten  using  the  definition  of  chemical  potential  a,,  for  species  /,  where  R„ 
is  the  universal  gas  constant. 


=G,”  +  /!,rto(a,), 

+  Ef=i  =  0  (/  =  l,2,...,iV).  (9) 

The  standard  Gibbs-energy  change  of  formation  for  species  i  is  denoted  which  is  equal  to  zero 
for  elements  in  their  standard  states.  The  activity  for  species  i  in  solution  is  given  by  cc,- ,  defined  in 
terms  of  the  equilibrium  constant  K  as 


where  the  activities  of  the  components  are  raised  to  flie  corresponding  stoichiometric  coefficients,  V(. 
For  an  ideal  gas  mixture  (.X;  <t) ,  -  1),  where  (|) ,  is  the  void  fraction. 


=  y;  =  x.^.p  =  p, 


fi  is  the  fugacity  and  X,  is  the  mole  firaction  for  the  ith  species.  For  liquids  and  solid  phases. 

In  (a,.)  =  In  (1  -  1/p), 

which  is  approximately  zero  for  large  pressure;  therefore,  6,  -  G®. 

There  are  N  equilibrium  equations  (equation  9),  one  for  each  species,  and  there  are  M 
material-balance  equations  (equation  8),  one  for  each  element,  a  total  of  AT  +  M  equations.  The 
unknowns  in  these  equations  are  the  n,’s,  of  which  there  are  N,  and  the  X^s,  of  which  there  are  M, 
a  total  of  + M.  Thus,  the  number  of  equations  is  sufficiait  for  the  determination  of  all  unknowns. 
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Numerical  experiments  were  performed  on  this  portion  of  the  code  as  well  to  independently 
verify  reaction  calculations  results  from  this  routine  with  known  results.  Principally,  known  gas 
phase  systems  wctb  investigated.  The  results  of  which  compared  exactly  with  those  of  the  NAS  A- 
Lewis  equilibrium  code  [25]. 

6.  Surface  Erosion  Results 

Due  to  current  development  programs,  the  first  application  of  the  coupled  ablation-chemical 
surface  reaction  model  has  been  directed  toward  direct-fire  systems.  The  model  is  used  to  compare 
the  relative  erosion  performance  in  the  M256  cannon  firing  die  M829A1  KE  projectile  and  an 
advanced  KE  projectile  concept  round.  Both  of  these  rounds  are  built  firom  similar  propellants, 
aldiough  the  flame  temperature  of  the  M829A1  propellant  is  lower  than  that  of  the  advanced  concept. 
The  thermal  performance  of  the  M829A1  KE  round  has  been  investigated  previously  [26]  with 
excellent  agreement  between  the  model  and  existing  experimental  data,  as  was  the  case  with 
numerous  previous  applications  [12, 13]. 

The  ARL  erosion  code  calculations  require  a  core  flow  thermochemical  ouqiut  fixrni  IBBLAKE, 
and  an  XKTC  output  file  including  the  spatially  varying  gas  pressure,  tonperature,  and  velocity  [27]. 
For  bodi  calculations,  the  barrel  was  assumed  not  to  have  a  chrome  layer.  Once  degraded,  the 
chrome’s  effectiveness  is  small  [28].  It  will  be  seen  that  this  assumption  possibly  has  some  validity 
for  the  advanced  KE  round,  while  for  the  less  viscerative  M829A1  round  diis  assumption  appears 
to  be  somewhat  excessive. 

Results  of  the  calculations  are  shown  in  Figures  5—10.  Figures  5  and  6  ^ow  surface  temperature 
predictions  for  bofli  the  M829A1  and  the  advanced  round  at  flnee  represratative  axial  locations, 
640  m,  1,050  mm,  and  1^50  mm  fiom  the  rear  face  of  the  tube  (RFT).  The  640-mm  location  is  near 
the  forcing  cone  region.  One  can  see  from  Figure  5  that  the  surface  just  reaches  the  specified  melt 
temperature  at  the  640-mm  axial  location,  while  at  the  more  distant  locations  down-bore,  the  surface 
does  not  reach  the  melt  temperature.  Experimentally,  bore  erosion  has  been  observed  at  both  the 
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640-mm  and  the  1,050-mm  locations,  but  not  at  the  1,350-mm  locations  in  tubes  having  fired  the 
M829A1  round  [28].  Looking  at  Figure  6,  what  is  first  noticed  is  that  the  melting  in  this  case  is 
much  more  severe  and  for  a  much  longer  duration  at  all  three  axial  locations  predicted.  Also  noted 
is  that  the  peak  temperature  for  this  round  is  reached  somewhat  later  and  longer  in  duration  than  fliat 
of  the  M829A1 .  This  is  due  to  the  different  propellant  grain  geometry,  which  was  designed  to  create 
just  such  an  effect  in  the  pressure  profile,  in  order  to  optimize  the  ballistic  performance  with  the 
given  constraints. 

The  erosion  for  the  M829A1  round  is  typically  located  near  the  forcing  cone  at  640  mm. 
However,  erosion  has  been  observed  initiating  somewhat  down-bore  at  the  1,050-mm  location  in 
some  guns.  This  “abnormal”  erosion  pattern  is  under  investigation  and  is  believed  to  be  due 
somewhat  to  the  interaction  of  the  forward  bore  rider  with  the  tube,  in  a  scraping,  or  perhaps  the 
aluminum  bore  rider  is  reacting  with  existing  atmospheric  oxygen  in  the  tube  forward  of  the 
projectile  [28].  This  would  increase  the  rate  of  chrome  removal  at  these  specific  locations  and  thus 
expose  the  underlying  bare  steel. 

Figures  7  and  8  show  the  predicted  resultant  surface  heat  fluxes  due  to  convection  and 
thermochemical  reactions.  Notice  that  in  both  cases  the  level  of  thermochemical  energy  imparted 
to  the  surface  is  on  the  order  of  1  %  of  the  total  energy.  This  is  intuitively  correct  because  convective 
heat  transfer  calculations  gave  very  good  agreement  with  experimental  data  [26].  The  M829A1  flux 
rises  and  decays  much  more  rapidly  than  that  of  the  advanced  round,  again  due  to  the  charge  design. 
The  advanced  round  shows  a  spike  at  the  peak  of  the  flux  that  is  believed  to  be  caused  by  the 
“slivering”  (the  point  at  which  the  grain  bums  through  its  web)  of  the  grain. 

Predicted  erosion  rates  are  shown  in  Figure  9  for  both  rounds.  The  erosion  rate  for  the  M829A1 
is  almost  an  order  of  magnitude  smaller  than  that  of  the  advanced  round.  Even  though  the  surface 
temperature  is  constant  during  melting,  the  erosion  rate  follows  the  heat  flux  contour,  including  the 
spike  at  the  peak  heating  rate.  The  erosion  rate  decreases  for  the  advanced  round  as  the  axial 
location  increases  from  the  RFT. 
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Figure  5.  Predicted  Surfoce  Temperatures  for  the  M829A1  KE  Round  in  a  Nonchromed 
M256  Cannon. 


Figure  6.  Predicted  Surface  Temperatures  for  an  Advanced  KE  Round  in  a  Nonchromed 
M256  Cannon. 
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Figure  7.  Predicted  Eneigy  Flux  Distribution  at  0.64  m  From  RFT  for  an  M829A1 KE  Round 
in  a  Nonchromed  M256  Cannon. 


Time  (ms) 


Figure  8.  Predicted  Energy  Flux  Distribution  at  0.64  m  From  RFT  for  an  Advanced  KE 
Round  in  a  Nonchromed  M256  Cannon. 
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Time  (ms) 

Figure  9.  Predicted  Erosion  Rates  for  an  Advanced  KE  Penetrator  at  Three  Axial  Locations 
and  the  Erosion  Rate  for  an  M829A1  KE  Penetrator  at  One  Axial  Location. 


Average  Erosion  for 
Advanced  KE/rnd  at  640mm 


^  Advanced  EE  at  640mm 
^  Advanced  KE  at  1060mm 
Advanced  KE  at  1860mm 
-  ‘  M829al  at  640mm 


Average  Erosion  for 
M829Al/rnd  at  640mm 


RSEEiaBSSi 


Time  (ms) 


Figure  10.  Predicted  Erosion  Depths  per  Round  for  an  Advanced  KE  Penetrator  at  Three 
Locations  and  tiie  Erosion  Depth  per  Round  for  an  M829A1  KE  Penetrator  at  One 
Axial  Location.  Also  Presented  Is  Averaged  Experimental  Erosion  Depths  per 
Round  for  Both  Penetrators. 
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Figure  10  reveals  the  amount  of  material  removed  per  round  fired,  for  botii  rounds.  The  M829A1 
erodes  much  less  at  the  forcing  cone  region  than  the  advanc^  round.  The  predicted  erosion  is  about 
16%  below  the  experimental  data  (average  of  39  rounds)  for  the  advanced  round  and  a  bit  high  for 
the  average  M829A1  round.  For  the  “abnormal”  erosion  described  previously,  the  M829A1 
calculations  were  within  the  data  range.  The  use  of  constant  physical  properties  in  the  calculations 
probably  accounted  for  some  of  this  discrepancy  [28, 29].  It  is  also  possible  that  for  the  M829A1 
rounds  more  chrome  existed  at  the  surface  near  the  forcing  cone,  thus  inhibiting  the  erosion. 


7.  Surface  Chemistry  Results 


Table  1  is  representative  of  the  chemistry  occurring  at  the  surface,  at  a  specific  location,  during 
the  peak  surface  temperature  time  step.  Both  the  unreacted  (input)  and  reacted  (output)  species  from 
the  thermochemical  calculation  are  presented.  A  total  of  34  species  were  considered  for  both  the 
M829A1  and  advanced  KE  penetrator,  although  many  of  the  potential  product  species  were  not 
formed.  A  user  may  choose  any  material  for  the  surface  or  any  additional  gas  products  fi:om  the 
propellant  combustion.  For  these  calculations  Oj,  COj,  Hj,  CO,  HjO,  Nj,  SOj,  and  NO  represented 
about  98%  of  the  combustion  products.  The  surface  was  considered  to  be  steel  consisting  of  Fe(A),* 
CrjOj,  NI(B),  and  C(GR).^  The  following  species  were  included  in  all  of  the  thermochemical 
calculations;  however,  they  converged  to  a  value  within  an  internally  specified  tolerance  on  the  mass 
fiaction  for  zero:  NI,  NIO,  NOj,  HNO,  HCO  RAD,  CRO,  CROj,  SO,  COS,  Sj ,  FEO(S),s  FEO(L),** 
FE,  FEO,  FEOjHj,  FE(L),  FE304(S),  OH,  OH+,  and  S. 

The  preliminary  equilibrium  chemistry  results  from  both  rounds  were  basically  the  same  because 
of  the  similar  propellant  formulation.  Due  to  the  chemical  equilibrium  assumption  and  the  affinity 
for  COj  formation,  carbon  and  oxygen  are  formed  from  the  “apparent  breakdown”  of  the  COj 
reactant  at  the  surface.  The  molecular  oxygen  reacts  with  the  molecular  hydrogen  to  make  water, 
while  the  carbon  forms  graphite.  The  added  FE(A)  goes  through  the  phase  transition  to  FE(C),*  and 
the  additional  NI(B)  melts.  During  this  time  step  and  at  this  location,  the  surface  is  regressing  and 
is  being  replenished  with  fiesh  “steel.”  The  CRA  in  the  steel  acts  inertly,  as  does  the  Nj  from  the 


*  A  -  Base-centered  cubic;  ^  GR  -  Carbon  graphite;  §  S  -  Solid;  **  L  -  Liquid;  *  C  -  Face-centered  cubic. 
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Table  1.  Representative  Chemical  Reactants  and  Products  During  Melting  of  the  Surface 


Energy  Released  (J/kg)  =  12218.7 

Time  (ms)  =  5.349 

Unreacted  Species  at  the 
Surface  (Input) 

Reaction  Products  at  the 
Surface  (Output) 

Species 

Mass  (kg) 

Fraction 

Mass  (kg) 

Fraction 

H2 

0.4553E-06 

0.1147E-01 

0.8000E-20 

0.2016E-15 

O2 

0.9061E-13 

0.2283E-08 

0.1270E-18 

0.3200E-14 

CO 

0.7709E-05 

0.1943E+00 

0.9855E-05 

0.2483E+00 

0 

u 

0.6656E-05 

0.1677E+00 

0.1747E-18 

0.4401E-14 

H20 

0.3950E-05 

0.9953E-01 

0.8019E-05 

0.2021E+00 

S02 

O.OOOOE4OO 

O.OOOOE+00 

0.2542E-18 

0.6406E-14 

N2 

0.1073E-05 

0.2704E-01 

0.1073E-05 

0.2704E-01 

NO 

0.2222E-11 

0.5599E-07 

0.1191E-18 

0.3001E-14 

C(GR) 

0.8020E-07 

0.2021E-02 

0.9766E-06 

0.2461E-01 

FE(A) 

0.1759E-07 

0.4432E-03 

0.2216E-18 

0.5585E-14 

NI(B) 

0.3299E-09 

0.8313E-05 

0.2330E-18 

0.5870E-14 

CR203(S) 

0.1634E-06 

0.4117E-02 

0.1634E-06 

0.4117E-02 

FE(C) 

0.1958E-04 

0.4934E+00 

0.1960E-04 

0.4939E+00 

NI(L) 

O.OOOOE+00 

O.OOOOE+00 

0.3299E-09 

0.8313E-05 

’>  _  2,;:. ' 

/-  .tv-'"'-'  .V  .  j 

Total  Mass  (kg) 
0.3969E-04 

Total  Mass  (kg) 
0.3969E-04 

Before  Sum  of  Mass 
Fractions  =  O.lOOOE+01 

After  Sum  of  Mass 
Fractions  =  O.lOOOE+01 
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gas.  One  interesting  point  is  that  the  graphite  formed  may  create  a  physical  property  modification 
if  it  dissolves  or  diffuses  into  the  surface  of  the  steel. 


8.  Discussion 

A  numerical  method  for  estimating  the  erosion  that  occurs  within  gun  tubes  due  to  high 
convective  and  chemicaUy  reactive  in-bore  heating  has  been  developed.  The  method  has  been 
successfully  applied  to  estimate  the  erosion  that  may  occur  during  the  firing  of  M829A1  and 
advanced  kinetic  energy  projectiles  fired  from  an  M256  cannon  that  has  the  chrome  chipped  away. 
While  chrome  is  present,  the  surface  does  not  reach  a  high  enough  temperature  to  enable  the  chrome 
nor  the  substrate  to  melt.  If  the  surface  is  gun  steel  (i.e.,  the  chrome  has  chipped  off),  then  the 
predicted  levels  of  erosion  are  roughly  in  agreement  with  experimental  observations  at  the  forcing 
cone.  The  prediction  for  the  M829A1  is  probably  high  due  to  the  absence  of  chrome  in  the  model, 
while  the  computed  values  for  the  advanced  KE  round  were  a  bit  low,  possibly  reinforcing  the 
assumption  that  the  chrome  was  not  very  effective  in  protecting  the  surface.  With  the  addition  of 
currently  excluded  physics,  such  as  surface  layers,  surface  defects,  gas  intrusion,  species  diffusion, 
and  in-wall  altered  material  properties,  the  predictions  should  improve. 


The  current  model  provides  the  U.S.  Army  with  the  ability  to  predict  the  effects  of  modifications 
to  the  propelling  charge  and  to  the  chemical  propulsion  technique  (i.e.,  liquid  or  solid  propellants). 
It  is  evident  that  small  chemical  formulation  changes  greatly  alter  the  erosion  behavior  [30]. 
Tracking  these  chemical  modifications  and  their  impact  on  the  erosion  behavior  is  now  possible. 
This  tool  may  be  especially  useful  in  estimating  the  erosive  behavior  before  the  development  of  new 
charges.  It  is  intended  to  be  used  as  a  platform  to  build  upon  to  include  the  mechanical  portion  of 
the  erosion  and  to  investigate  othar  physics  of  interest  such  as  the  effect  of  wear-reducing  additives 
and  surface  treatments.  Potentially,  newer  and  more  direct  mitigation  techniques  may  be  identified 
once  the  causality  of  a  specific  erosion  problem  is  determined.  Also,  this  model  has  been  written 
to  incorporate  boundary  conditions  from  future  interior  ballistic  codes  such  as  NGEN  [31].  NGEN 
will  be  fully  turbulent,  enabling  the  turbulent  reacting  boundary  layer  to  be  computed.  Tliis  coupling 
will  strengthen  the  mass  and  heat  transfer  physics  and  enable  more  confidence  in  the  solutions. 
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